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Abstract 

The status of lattice QCD investigations at high temperature is re- 
viewed. After a short introduction into thermal QCD on the lattice 
we report on the present understanding of the phase diagram and the 
equation of state, in particular in presence of dynamical quarks. We 
continue with a discussion of various screening lengths in the plasma 
phase including results from dimensionally reduced QCD. This is fol- 
lowed by summarizing lattice data on quark number susceptibilities 
and spectral densities, both of which are of immediate relevance to the 
interpretation of heavy ion experiments. A major section is devoted 
to presenting simulations of QCD at small yet phenomenologically 
important values for the baryon density. 
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1 Introduction 



Quantumchromodynamics (QCD) is the theory describing the strong inter- 
actions, carried by gluons, which confine quarks into hadrons. It is the 
fundamental theory underlying all nuclear physics, and as such responsible 
for many salient features of matter as we know it today. A key property of 
the theory is asymptotic freedom, according to which the coupling strength 
decreases with the energy transfer of an interaction. One consequence is that 
at high enough energies of a few GeV the theory is tractable by perturba- 
tion theory, which leads to its experimental verification in e.g. deep inelastic 
scattering. On the other hand, on the low energy scales of hadron physics 
the coupling is strong and perturbation theory invalid. The only known non- 
perturbative and first principles method to compute QCD predictions in this 
regime is by simulations of lattice gauge theory, whose results are beginning 
to give a quantitative description of the hadron spectrum. 

At a critical temperature T c ~ 200 MeV, QCD predicts a phase transition 
to a deconfined plasma of quarks and gluons. At the same temperature the 
weakly broken chiral symmetry, responsible for the existence of light pions, 
gets restored. Such temperatures were realized in the early universe, which 
during its cooling expansion passed through the plasma phase and the quark 
hadron transition on its way to its present state. We are currently witnessing 
exciting attempts to recreate this primordial plasma and its subsequent cool- 
ing through the phase transition in heavy ion collision experiments at RHIC 
(BNL) as well as SPS and LHC (CERN). These studies will have a bearing 
far beyond QCD in the context of early universe physics. Many prominent 
features of the observable universe, such as the baryon asymmetry or the 
seeding for structure formation, have been determined primordially in hot 
plasmas described by non-abelian gauge theories. The QCD plasma is pro- 
totypical for those as it is the only one we can hope to produce in laboratory 
experiments. Moreover, for high densities and low temperatures more exotic 
non-hadronic phases like a color superconductor have been predicted, and 
there is a chance for such physics to be realized in the cores of compact stars. 

For such applications we need to know how the properties of QCD change 
under extreme conditions of temperature and density. This entails a deter- 
mination of the QCD phase diagram and the associated critical properties, a 
quantitative understanding of the equation of state, the way in which hadron 
properties get modified as well as the nature of the lightest excitations in 
non-hadronic phases. Answers to those questions may also help to elucidate 
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mechanisms of confinement and chiral symmetry breaking. 

In this review we collect data from lattice simulations of thermal QCD and 
try to place them into their physics context. While avoiding technical issues 
as far as possible, our aim is to convey a picture of what we have learned 
about the QCD transition and plasma, even though in places this picture 
is still incomplete. When the critical temperature is crossed, light degrees 
of freedom appear to be released, in accordance with the naive expectation 
of deconfmement. However, the resulting plasma is not compatible with the 
perturbative picture of a free parton gas for any temperatures of interest. 
Rather, the system displays strong residual intercations caused by soft modes 
in the gauge sector. 

Static equilibrium physics is well under control in the continuum limit for 
the pure gauge theory without quarks, whereas calculations with dynamical 
quarks are limited to m n > 300 MeV, and even those are still affected by 
cut-off effects. Nevertheless, these systematic errors can be reduced step by 
step in additional calculations to come. 

Less advanced and not yet fully quantitative are the calculations in two 
other areas of interest. Since lattice QCD works in euclidean space-time, it is 
much harder to obtain a real time description of dynamical processes. This 
field is in its infancy and we will report on a first step taken in that direction. 
In the last two years significant progress has been made towards simulating 
the regime of small baryon densities relevant for heavy ion collisions, which 
was entirely impossible before then. While this also is a rather young field 
of research with a long way to go, there now are first results for the (fi, T) 
phase diagram and the search for a potential critical endpoint of the decon- 
fmement transition. Before presenting lattice results available today, we give 
a brief introduction to the theoretical set-up of lattice thermodynamics and 
a discussion of its current limitations. 

2 QCD Thermodynamics on the Lattice 

Detailed introductions to thermal field theory and lattice field theory can 
be found in pQ and j2], respectively. The QCD grand canonical partition 
function is given by 

Z(V,fi,T) = Tr ( e -^-^)/ T ) = J DAD^jD^j e - s ^e~ s ^ A ^ , (1) 
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with the euclidean gauge and fermion actions 

Its thermodynamic parameters are the temperature T, the volume V and 
the chemical potential /i for quark number Q (the one for baryon number is 
fiB — 3/x), while the QCD action depends on the number of quark flavors Nf, 
their masses and the gauge coupling g. In the following we will consider 
mostly two and three flavors, and always take m u = rrid- The case m s = m u ^ 
is then denoted by Nf = 3, while Nf = 2 + 1 implies m s ^ m u ^ d . 



l/T 

= J dx J d 3 x -Tr F^F^, 
o v 

= J dx J d 3 x J2i! f {y^D^ + ml - fi^tpf. (2) 



2.1 The Lattice Formulation 

The theory is discretized by introducing a euclidean space-time lattice L 3 x iV t 
with lattice spacing a, such that volume and temperature are 

The fermion fields live on the lattice sites x, whereas the gauge fields are 
represented by link variables U^(x) G SU(3) connecting the sites. After a 
suitable discretization of the actions, the Gaussian integral over the quark 
fields can be performed leading to the lattice partition function 

Z{L, an, N t ; (3, N f , am[) = f DU J] det M{n)e~ s '^, (4) 

/ 

where M(fi) = 7 M (-D M + m/ — //70) is the fermion matrix and S g [U] denotes 
the gauge action with lattice gauge coupling (3 = 6/g 2 . For \i = 0, the theory 
in this form is amenable to a stochastic calculation of expectation values by 
Monte Carlo methods. Simulations at \x 7^ have only begun to become 
possible in certain limits, as will be discussed in detail in Section El 

The lattice spacing depends on the gauge coupling, a = a(g), and is set 
by some physical (zero temperature) quantity like a hadron mass. The phys- 
ical temperature is then obtained via Eq. © as T/m# = l/(arriHN t ). In 
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general, calculations are affected by finite size and cut-off effects, the latter 
depending on the particular discretization chosen. The task then is to com- 
pute observables for various volumes and lattice spacings, and extrapolate 
to the thermodynamic (V — ► oo) and continuum limits (a — ► 0, Nt — > oo), 
while keeping T and physical parameters fixed. 

In practice calculations are severely limited by computational resources, 
the most expensive part being the evaluation of the fermion determinant. 
For this reason it has often been set to 1 in the past, corresponding to the 
quenched approximation which neglects all quark loop contributions. In dy- 
namical simulations at finite temperature, the Wilson and staggered fermion 
actions are used predominantly, as they are the cheapest to simulate. An- 
other choice with better chiral properties are domain wall fermions, while 
dynamical overlap fermions so far have not been employed. 

Typical lattice sizes of current dynamical simulations are 16 3 x 8 or 32 3 x 4, 
but for many exploratory studies smaller lattices have to make do. In order 
to have reasonably small finite size and cut-off effects, the correlation length 
£ of a typical hadronic state has to fit comfortably in the box while being 
much larger than the lattice spacing, 

a < £ < aL. (5) 

Around the deconfinement transition T c ~ 200 MeV ~ (lfm) -1 , the temporal 
size N t = 4 — 8 implies lattice spacings of a ~ 0.1 — 0.3 fm, allowing for box 
sizes of aL ~ 1.5 — 3 fm. In the chirally broken phase the right inequality then 
constrains the lightest feasible quark masses, implying that the physical pion 
cannot be accomodated on such lattices. Typical calculations correspond to 
m„ ~ 300 MeV. 

Above T c meson correlation lengths typically scale as £ ~ 1/T so that the 
constraint of Eq. (jSJ) reads 

w t <<l<< w t - (6) 

Cut-off effects are then controlled by 1/Nt independent of temperature. These 
cut-off effects can be rather severe for certain thermodynamic quantities. The 
reason is that thermal distribution functions are peaked for modes with mo- 
menta p ~ T of the order of the cut-off, or in lattice units ap ~ 1/Nt. Control 
over these effects thus warrants the use of improved actions, which are aug- 
mented by additional terms irrelevant in the continuum limit and tuned such 
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that, at finite lattice spacing, they subtract cut-off effects to a given order 
(a n ) in the lattice spacing. In addition, cut-off effects can be smaller on 
anisotropic lattices, for which a finer spacing in the temporal direction is 
chosen, at < a s , corresponding to a larger number of points N t . However, 
soft gluonic modes scale as £ ~ l/g 2 (T)T, upon which l/iV t — > g 2 (T)/N t in 
Eq. (JHJ), thus requiring large aspect ratios L/N t in both cases, which limits 
simulations of very large temperatures. 

2.2 Effective High T Theory: Dimensional Reduction 

At large T, when the gauge coupling g(T) is sufficiently small, a hierarchy 
between different relevant scales of thermal QCD develops, 

2ttT > gT > g 2 T. (7) 

The lowest non-vanishing bosonic Matsubara mode ~ 2ttT is characteristic 
for non- interacting particles. The dynamics generates the Debye scale mg~ 
gT, over which color electric fields are screened, and its non-perturbative 
analogue tum ~ g 2 T for color magnetic fields [3]. 

For physics on scales larger than the inverse temperature, |x| ~ 1/pT ^> 
1/T, this allows an effective theory approach in which the calculations are 
factorized: integration over the hard modes may be performed perturbatively 
by expanding in powers of the ratio of scales gT/(27iT) ~ g/(2ir). This 
includes all non- zero Matsubara modes, in particular the fermions. It results 
in a 3d effective theory for modes ~ gT and softer, 

S eff = J d 3 x {^Tr(iV^) + Tr(A^o) 2 + m%Tr(A 2 ) + A 3 (Tr(^) 2 | , (8) 

and is known as dimensional reduction jjj. Since 4d euclidean time has 
been integrated over, A now appears as an adjoint scalar, and the effective 
parameters are functions of the original ones, g\ = g 2 T, rriE{N, Nf, g, mQ ~ 
gT,X 3 (N,N f ,g,mf) ~ g A T. 

Discretization and simulation of this reduced problem is easy. Without 
fermions and one dimension less, much larger volumes and finer lattices can 
be considered. Moreover, 3d gauge theories are superrenormalizable and the 
coupling scales linearly with lattice spacing. Hence, very accurate continuum 
limits can be obtained and systematic errors from simulations are practically 
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eliminated. However, the reduction step entails two approximations: the per- 
turbative computation is limited to a finite order in g and neglects higher- 
dimensional operators, which are suppressed by powers of the scale ratio. 
The reduction step has been performed up to two-loop order [5] at which 
parameters have relative accuracy 0(g 4 ), while for correlation functions the 
error is jB] SC/C ~ 0(g 3 ). In the treatment of the electroweak phase transi- 
tion, this is less than 5% [7], for QCD it depends on the size of the coupling 
g{T). As we will discuss in Section it can be non-perturbatively checked 
that the effective theory accurately describes static correlation lengths at 
temperatures as low as T > 2T C , thus allowing for a straightforward treat- 
ment of very large temperatures as well as detailed dynamical investigations 
in the plasma phase. 

3 The (m Uj d, m s ) Phase Diagram 

On quite general grounds it is expected that at high temperature or baryon 
density QCD undergoes a transition from hadronic matter to the quark- 
gluon plasma. However, whether the transition is characterized by a truly 
singular behavior of the partition function leading to a first or second order 
phase transition, or whether it merely is a crossover with rapid changes in 
some observables, crucially depends on the values of the quark masses. At a 
transition the details of the microscopic interactions are negligible compared 
to the dominant long range correlations, which are determined by global 
symmetries. Expectations on the nature of the transition have thus been 
derived from studies of simpler systems with the same global symmetries as 
QCD |8], see Figure HJ 

The method to locate a transition usually is to look for rapid changes 
of order parameters like the Polykov loop L (in the pure gauge case) or the 
chiral condensate (in the chiral limit) and for peaks of their fluctuations, 
i.e. susceptibilities as e.g. xl = V((L — (L)) 2 ). The locations of these 
changes or peaks define (pseudo-) critical couplings /3 C , which can be turned 
into temperatures through the knowledge of a zero temperature quantity like 
a hadron mass am#(/3 c ) at the same coupling, T c /m# = l/(N t amH(Pc)- The 
nature and critical properties of the transition can be obtained from scaling 
dependencies of various quantities on volume or external fields like the quark 
mass. 

In the 577(3) pure gauge theory, the expected first order transition was 
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Figure 1: The phase diagram, expected (left) and lattice data (right), in the 
plane of strange and degenerate u,d quark masses. 



verified numerically |9 3 quite a while ago. In the chiral limit of two-flavor 
QCD the relevant global symmetry group most likely is SUr(2) x SUl(2) ~ 
0(4). In analogy to the 0(4) a model, one therefore expects [TU] a second 
order phase transition with 0(4) critical properties in the chiral limit. At 
non-vanishing quark masses the chiral symmetry is broken explicitely and 
a crossover is expected. All existing lattice calculations jTTJ [T2J H31 UH EES] 
support a continuous transition. However, 0(4) values for the critical expo- 
nents could not yet definitively be identified. Alternatively, the anomalous 
XJ a (1) could become restored effectively through the disappearance of topo- 
logically non-trivial zero modes of the Dirac matrix at or even below the 
critical temperature rendering the transition first order JU] ■ This case would 
be reflected in a degeneracy ^B] between the pion and the scalar isovector 
meson (<5/a ), which does not seem to be realized jT7j. 

For three degenerate massless quarks, again from symmetry reasons, one 
expects a first order transition. The first order property of the transition 
should extend from the chiral limit to an end point at some non-zero critical 
quark mass m c , where the transition becomes continuous and is expected to 
be in the Ising universality class [TH]. A recent lattice simulation [TH] supports 
this picture even with respect to the universality class. The location of the 
end point, however, is not a universal property and the current estimates 
[T9"l 12*0*1 |2"T] vary between 190 and 260 MeV for the corresponding pseudoscalar 
mass, depending on the particular fermion action chosen. 
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For non-degenerate quarks, i.e. a strange quark heavier than the two light 
ones, m s ^ m u ^ d , the end point becomes an end line separating the first order 
from the crossover regime, see Figure H in the vicinity of the degenerate 
three-flavor end point the slope of the line can be obtained from a Taylor 
expansion, m c s = m c + 2(m c — m U) d). In the currently accessible quark mass 
range, this relation is consistent with all available lattice data [!$] I2U1 122] 123] . 
see Figure ^ If one assumes that such a linear extrapolation holds down to 
the physical light quark mass values, the resulting ratio m c s /m^ s ~ 5 — 10 
would put the point of physically realized QCD, m p s hys / 'vnP^Y ~ 20 (see e.g. 
for current lattice results), into the crossover regime of the phase diagram. 

The location of the phase transition, i.e. its critical temperature, has been 
known rather precisely for quite a while in the pure SU (3) gauge theory [23 
I2E] . It has since then been confirmed by studies which explored the efficiency 
of various improved actions [2E1 1213 I2Z] , or anisotropic lattices [30] |3T] , in 
reducing finite lattice spacing effects. The number is most readily given in 
terms of the dimensionless ratio to the string tension, y/a ~ 425MeV, as 



which was obtained as a weighted average over all available lattice data. 

For QCD including light dynamical quarks the current results for the crit- 
ical temperature are summarized in Figure |2J The left panel contains data 
(221 1221 1211 HUB inn] for two flavors, all obtained with improved actions. Dis- 
crepancies between Wilson and staggered quarks seen in earlier computations 
with standard discretizations are greatly reduced. T c is plotted versus the 
zero-temperature pseudoscalar mass as a measure of the quark mass, both in 
units of the vector meson mass. The strong decrease of the ratio T c /my with 
increasing quark mass is due to the rising vector meson mass which diverges 
in the infinite quark mass limit, mps/my — > 1. Close to the chiral limit 
two conflicting quark mass dependencies should influence the ratio T c jmy\ 
the vector mass should increase linearly with the quark mass, leading to a 
decrease of the ratio with y/mps- This is counteracted by the quark mass de- 
pendence of T c itself, rising proportional to m^ 136 , which amounts to ~ rnp^s 
for 0(4) values for the critical exponents. Thus, at small enough quark and 
pseudoscalar mass, respectively, a linear decrease towards the chiral limit is 
expected to win. This effect is not yet seen. A scale for T c much less affected 
by the quark mass is the string tension which is used on the right panel of 
Figure |2] Indeed, this figure shows that the critical temperature is decreasing 




(9) 
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with decreasing quark mass. The decrease appears to be compatible with a 
linear dependence on the pseudoscalar mass as expected for Nf = 2. The 
linear dependence, however, also seems to hold for three degenerate quarks. 
Note also, that T c is reduced half way down from the quenched value, shown 
as the band to the very right of the figure, already at a pseudoscalar mass 
as large as 1500 MeV. It is thus not the chiral dynamics which controls the 
critical temperature, and one may speculate that a resonance gas picture is 
more appropriate to describe the thermodynamics close to T c . The overall 
quark mass dependence for both two and three flavors is rather weak. The 
difference between Nf = 2,3 is small and amounts to about 20 MeV, inde- 
pendent of the quark mass. Extrapolating the results to the chiral limit, one 
obtains 

(171 ± 4) MeV, clover-improved Wilson 

fermions [33] 
(173 ± 8) MeV, improved staggered 

fermions 32] 
(154 ± 8) MeV, improved staggered 
fermions [32] 

where m p has been used to set the scale. Although the agreement between 
staggered and Wilson quarks is striking, one should keep in mind that the 
errors are statistical only and do not account for the as yet unknown sys- 
tematic errors originating from a non-vanishing lattice spacing. One might 
hope, though, that these are small since improved actions have been used in 
both simulations. 

4 Equation of State 

Energy density e(T) and pressure p(T) as a function of temperature are cer- 
tainly among the most fundamental thermodynamic quantities of the quark 
gluon plasma governing, e.g., the time evolution of the plasma once being 
created in heavy ion collisions. At high temperature the relevant partonic de- 
grees of freedom have momenta of order irT which are strongly affected by the 
UV cut-off introduced through the finite lattice spacing. It is thus important 
to gain control over the discretization effects and carry out the continuum 
limit a — > 0. Again, this was first accomplished in the pure gluon theory 
|25j . later confirmed in studies utilizing improved actions and/or anisotropic 



2 - flavor QCD : T c = 

3 - flavor QCD : T c = 



11 




Figure 2: Left: T c in units of the vector meson mass for Nf = 2, for a variety 
of improved actions [H2J EH EH ESI EE] • Right: Nf = 2, 3 for improved (p4) 
staggered quarks, compared to Nf = 2 standard (std) simulations [32] • 



lattices j2*Sl 12*^ I3T] . As calculations in the ideal gas limit of lattice perturba- 
tion theory reveal, the discretization effect is even more pronounced in the 
presence of dynamical quarks. Since on the other hand the signal for these 
quantities vanishes proportional to a 4 the use of improved actions, which are 
designed to reduce the UV cut-off effects, seems mandatory. Moreover, the 
simulations in the quenched approximation have shown that this improve- 
ment helps to extract the continuum limit also in the intermediate T range 
between 2 and 4 T c . 

The results of a computation of the pressure with an improved action in 
the gauge as well as in the fermion sector [37] are shown in Figure El The data 
have been obtained for Nf = 2, 3 with (bare) mass m q /T = 0.4 as well as for 
Nf = 2 + 1 with a heavier mass m s q /T = 1. For comparison the figure includes 
the continuum extrapolated quenched result. The figure shows a rapid rise 
of the pressure in a transition region. The critical temperature as well as the 
magnitude of p/T A reflect the number of degrees of freedom liberated at the 
transition. In fact, when the pressure is normalized to the appropriate Stefan- 
Boltzmann limits shown as arrows in the figure, the function p(T/T c )/psB 
turns out to be almost flavor independent. It is thus well described by 

^=(l6 + yiV / )^x / (r/r c ), (10) 

where the ideal gas limit only is modified by an apparently flavor indepen- 
dent function f(T/T c ). At this stage it should be pointed out, though, that 
the full QCD results have not yet been extrapolated to the continuum limit. 
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Figure 3: Left: flavor dependence of the pressure for N t = 4 lattices com- 
pared to a continuum extrapolated pure gauge result [37]. Right: energy 
density for Nf = 2 improved Wilson quarks on N t = 4 (filled symbols) 
and N t = 6 (open symbols) 39\. Marks on the right side denote Stefan- 
Boltzmann limits. 

However, the experience gathered perturbatively and in the quenched ap- 
proximation leads one to expect that the finite a effects will distort these 
findings by not more than 10%. 

The most recent results on the energy density are shown in Figure [3] 
to the right. Note that the dependence on the quark mass appears to be 
rather weak over a wide range of values. The data has been obtained for 
Nf = 2 improved Wilson quarks jSH] on Nt = 4, 6 lattices. However, in 
the particular discretization chosen the improvement does not reduce the 
infinite temperature finite lattice spacing effects. This is visible in the figure 
by the big difference in the Stefan-Boltzmann limits computed on the two 
lattices with different N t values. On the other hand, the cut-off effects are 
not important close to T c as a comparison between N t = 4 and 6 results 
reveals. This is expected since in this temperature regime the correlation 
length is large and infrared modes dominate. Thus, a comparison of the 
different fermion discretizations is meaningful. In fact, improved staggered 
data are available jlU] and lead to a consistent estimate for the critical energy 
density of e(T c ) ~ (6±2)T C 4 , which also agrees with the value originating from 
an earlier simulation using the standard staggered action |41j . 

In all cases, up to the highest temperature investigated, pressure and 
energy deviate substantially from the ideal gas behavior. The deviation is 
too big to be reproduced in ordinary high temperature perturbation theory 
which converges badly in this temperature range [12]. A more refined tech- 
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Figure 4: The pressure of pure gauge theory from dimensional reduction, 
with an as yet undetermined constant e$. From |4*4] . 

nique as hard thermal loop resummation [3B] gets closer to the lattice results. 
The question of the approach to the ideal gas limit at higher temperatures 
can also be addressed by dimensional reduction, which splits the thermody- 
namic functions into perturbatively calculable parts and a remainder to be 
simulated in the effective theory Eq. (JHJ) 0Bj- Unfortunately, this requires 
determination of an integration constant 6q by a four-loop calculation, which 
has not been completed yet. However, results (Hj treating e like a free pa- 
rameter and matching it to the 4d simulations are shown in Figure 0] If the 
eventual value for eo lies in the right range, the whole picture makes sense and 
explains how ideal gas values are only attained at asymptotic temperatures. 



5 Screening Masses 

Essentially all static equilibrium properties of a thermal quantum field theory 
are encoded in its equal time correlation functions. Even though these quan- 
tities may not be directly accessible in heavy ion collision experiments, their 
theoretical knowledge provides us with the relevant dynamical length scales 
in the plasma, from which conclusions about the acting degrees of freedom 
and their physical effects may be drawn. The connected spatial correlation 
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functions of gauge- invariant, local operators A(x), 

C(|x|) = (A(x)A(0)) c ~e- M W, (11) 

fall off exponentially with distance. The "screening masses" M are the eigen- 
values of the spacewise transfer matrix of the corresponding lattice field 
theory, and classified by its symmetries. Because of the shortening of the 
euclidean time direction at T > 0, the continuum rotation symmetry of the 
hypertorus orthogonal to the correlation direction is broken down from 0(3) 
to 0(2) x Z(2), and its appropriate subgroup for the lattice theory is D\. 
The irreducible representations and the classification of operators have been 
worked out for pure gauge theory [13 HH] as well as for staggered quarks [H] . 
Physically, the screening masses correspond to the inverse length scale over 
which the equilibrated medium is sensitive to the insertion of a static source 
carrying the quantum numbers of A. Beyond 1/M, the source is screened 
and the plasma appears undisturbed. Technically, the computation of these 
quantities is equivalent to spectrum calculations at zero temperature. 

5.1 Hadronic Screening Masses 

Figure (left) shows results for the lowest lying screening masses correspond- 
ing to glueball operators around T c . Comparison with dimensionally reduced 
results (cf. Figure \§§ shows that the pure gauge results are close to contin- 
uum physics. In the range 0.8T C < T < T c , the lowest scalar screening mass 
is observed to be roughly 20% lower than the lightest scalar glueball at zero 
temperature, M{T) / M G {T = 0) ~ 0.8. At T c a sharp dip is observed, after 
which the screening masses appear to be proportional to T. Screening states 
with larger masses show the same qualitative behavior above T c . 

Apart from operators A made from pure glue in Eq. (jllj) . also spa- 
tial correlations of meson operators have been investigated, both in the 
quenched approximation as with various numbers of dynamical fermions 
US EOJ Ell E2J E31 EU ESI ESI EZI The picture which has emerged so far 
is illustrated by some of the available data shown in Figure El (right). Below 
T c , the screening masses have shown neither a marked temperature depen- 
dence nor, correspondingly, a drastic difference to zero temperature masses. 
At temperatures above T c spatial (as well as temporal) correlation functions 
reflect the restoration of the chiral SU^Nf) x SUn(Nf) symmetry. In par- 
ticular, the vector and axial vector channel become degenerate independent 
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Figure 5: Left: Screening masses for the pure gauge theory, corresponding to 
the continuum 0+ + (circles) and 0l~ (squares) channels. From [IB]- Right: 
Mesonic screening masses in the quenched chiral limit. Below T c , M p /T c 
jSH E2] is plotted, the line denotes the T = value. Open Wilson symbols 
denote anisotropic lattices Staggered pion data are extrapolated to 

a = j56j, Wilson and staggered rho are from N t = 8 [HI] and 16 [S2] 
lattices. The free quark limit has not been corrected for finite volume effects. 



of the discretization and of the number of dynamical flavors being simulated. 
Moreover, the pion ceases to be a Goldstone boson and acquires a screening 
mass. A degeneracy is further observed within errors in the pseudoscalar 
and isoscalar scalar (cr// ) channel although the latter, for technical reasons, 
is difficult to access on the lattice. Nevertheless, screening masses obtained 
from fits to correlation functions [HZ] and susceptibilities (at finite lattice 
spacing) [58] lead to a consistent picture. At T c the anomalous Ua(X) most 
likely is not yet restored effectively [T7j, see also [SHI EH EE EH] > as has been 
mentioned already in Section El 

At high temperature, the screening masses are expected to approach the 
free quark propagation limit, M — > 2irT. In fact, already at temperatures of 
about 1.5 T c , the results for the vector channel are not far from this value. 
If the finite volume effects for free quark propagation are taken into account, 
the deviations amount to about 15 % and decrease only slowly with temper- 
ature. Moreover, at temperatures > 1.5T C , in the Wilson discretization pion 
and rho become nearly degenerate already at finite lattice spacing. Earlier 
staggered simulations had found discrepancies here. However, a recent paper 
[HS] reports that also in the staggered discretization a degeneracy of pseudo- 
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Figure 6: Left: Pure gauge screening masses at 2T C from simulations of the 
4d (circles) [46J and the dimensionally reduced [HU E2j theories, where filled 
squares correspond to 0j; + : Tr(F^ 2 ) : m glue , open squares represent 0| + : 
Ti(Al) : m Ao ,0z + : Ti(A F 12 ) : m D and 0±~ : Tr(Ajj). Right: Temperature 
dependence of the two lowest 0| + and the lowest 0z + states 153]. 



scalar and vector is reached in the continuum limit. Quenching effects are 
found to be below 5% for T > T c 60J. 

In the plasma phase, screening masses can also be investigated within 
the dimensionally reduced theory, cf. Section 12.21 In this framework they 
correspond to the spectrum of the transfer matrix for the 3d theory. The as- 
sociated Hamiltonian respects SO (2) planar rotations, two-dimensional par- 
ity P, v4o-reflections R, and the symmetry is again 50(2) x Z(2) x Z(2) = 
0(2) x Z(2). Remember however that in this setup one is interested in soft 
modes, while the Matsubara frequency ~ 2-rT represents the UV cut-off for 
the effective theory. 

Figure El displays numerical results for the lowest screening masses for the 
SU(2) and SU(3) pure gauge theories, obtained through simulations of both 
the 3d effective theories and the full 4d theories, in various quantum number 
channels. Even at a temperature as low as T ~ 2T C , good agreement is ob- 
served except in one case. Note that the error bars are only statistical, while 
both calculations have systematic errors: the effective theory has trunca- 
tion errors, while the 4d simulations are not infinite volume and continuum 
extrapolated. Estimating such effects at 10-15%, we conclude that within 
current accuracy dimensional reduction gives a quantitative description of 
screening lengths for temperatures T > 2T C . The same conclusion is reached 
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with other observables like the spatial string tension [61| . static potentials 
or gauge fixed propagators (cf. Section lo~3|) . For larger temperatures the 
M/T get lower and the spectrum denser, as shown in Figure El (right) for 
a few states. A large number of such screening masses in many quantum 
number channels for various numbers of flavors has been computed in [U2"| . 

5.2 How Perturbative is the Plasma? 

The approach of dimensional reduction is particularly valuable in disentan- 
gling contributions from different degrees of freedom, thanks to accurate 
mixing analyses, as well as for treating larger temperatures T ^> T c which 
cannot be reached by 4d lattices. This can be used to inspect to what ex- 
tent the plasma behaves perturbatively. Figure |B] then teaches us that for 
any reasonable temperature the largest correlation length of gauge-invariant 
operators belongs to the Aq degrees of freedom and not to the A4, in con- 
trast to the naive parametric ordering, Eq. ((7j). Only asymptotically is the 
perturbative ordering restored. 

Of particular phenomenological relevance for the QCD plasma is the De- 
bye mass, whose inverse gives the length scale over which colour-electric flux 
is screened. Its expansion in powers of coupling constants reads 

m D = m D ° + ^ In ^ + c N g*T + OijfT), (12) 

with the leading order perturbative result m D ° = (N/3 + Nf/Q) l l 2 gT. At 
next-to-leading order ~ g 2 T, only a logarithm can be extracted |66j . whereas 
the coefficient cat is entirely non-perturbative. 

A gauge invariant definition of the Debye mass has been suggested in |HZj , 
according to which it corresponds to the mass of the lightest gauge invariant 
screening state odd under euclidean time reflections. In the 3d effective 
theory the latter is replaced by the scalar reflection R, and the Debye mass 
thus corresponds to the 0z + ground state ~ TrA Fi 2 . 

The coefficient can be measured separately from the exponential decay 
of a Wilson line in a 3d pure gauge theory [HZ], which has been performed 
in |68j . This allows to disentangle the contributions from different scales, as 
shown in Tabled The 0(g 3 T) corrections are less than 30% even at temper- 
atures as low as T = 2T C , and they disappear entirely for asymptotically large 
temperatures. On the other hand, the dominant scale is again ~ g 2 T for all 
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Table 1: Contributions of the first two, the third and higher order terms in 
Eq. (TJJ) to the total m D £2|. 



SU(3), N f = 


m D /g 2 T 


(1 + 2) 


c 3 


0(g 3 T)/(giT) 


T = 2T C 


1.70 (5) 


0.514 


1.65(6) 


-0.46(6) 


T = 10 n T c 


3.82 (12) 


2.165 


1.65(6) 


0.00(12) 



temperatures of interest, in contrast to the naive expectation ~ gT. Only at 
asymptotically high temperatures is the perturbative picture restored. These 
findings are compatible with and explain the behavior of the pressure density, 
which attains its ideal gas limit only at asymptotic temperatures. There the 
effect is less pronounced, since the pressure is dominated by hard modes and 
soft modes provide corrections, whereas the screening masses discussed here 
couple exclusively to soft modes. 

5.3 Gluon Screening Masses 

With energy and pressure densities as well as the quark number suscepti- 
bilities rising rapidly across the transition, it seems clear that some light 
constituents are liberated in the deconfined phase. A longstanding ques- 
tion concerns the nature of these degrees of freedom and their properties. 
Moreover, spatial correlators of singlet operators as discussed in the previous 
subsections, and in particular a Debye mass defined by one of them, are not 
directly related to phenomena like J/ ^-suppression. These are caused by 
color charged intermediate states, i.e. some kind of constituent. In perturba- 
tion theory one therefore defines the Debye mass in analogy to QED as the 
pole in the A -propagator [66J, with a perturbative expansion as in Eq. (JHJ, 
which can be proved to be gauge invariant order by order On the other 
hand, at moderate temperatures T > T c interactions are still strong, and a 
purely perturbative parton picture is surely not appropriate. Similarly, the 
propagator of magnetic fields is expected to develop a magnetic mass on 
the scale ~ g 2 T, which however even at asymptotic temperatures is entirely 
non-perturbative [3]. 

Unfortunately, lattice studies of gluon propagators [70] are hampered by 
several problems. It is difficult to fix a gauge uniquely and avoid the prob- 
lem of Gribov copies [7T]. Moreover, most complete gauge fixings, like Lan- 
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Figure 7: The electric screening (Debye) mass in SU(2) gauge theory. The 
line represents a fit to 4d simulations [73] , the symbols show simulations of the 
dim. red. theory with different matching procedures. Triangles correspond 
to two-loop perturbative matching. From 



dau gauge, violate the positivity of the transfer matrix, thus obstructing a 
quantum mechanical interpretation of the results. However, recently it was 
shown analytically [72 J that gluon correlators can be made gauge invariant by 
dressing them with appropriate (gauge fixing) functionals. These correlators 
fall off exponentially with eigenvalues of the Kogut-Susskind Hamiltonian in 
the presence of sources [72], thus probing the energy of a gluon coupled to 
sources. The spectrum is gauge invariant, while matrix elements depend on 
the particular functionals chosen. This is true for all gauge fixings which are 
local in time (to preserve the transfer matrix), such as the Coulomb gauge. 

Numerically, in detail only SU(2) pure gauge theory has been investi- 
gated so far. Figure [7J shows results obtained for the electric screening mass 
obtained both in 4d simulations as well as in simulations of the dimension- 
ally reduced theory, where a variety of gauges was found to agree with the 
Coulomb gauge results. Like for the hadronic screening masses, with two-loop 
matching good agreement between full and effective theory for temperatures 
down to T « 2T C is observed. Equally consistent is the finding that the lead- 
ing ~ gT contribution to the electric mass is subdominant for temperatures 
up to O(10 4 T c ) [TJ. 

Extracting the magnetic mass appears to be more difficult technically. In 
Landau type gauges, where positivity is violated, severe finite volume effects 
are observed [731 UM- m Coulomb type gauges the spatial non-locality of 
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the fixing functionals leads to an exclusive projection on torelonic states 
[72]. However, by introducing explicit local fields to mimic static sources 
and extrapolating to their infinite mass limit, the magnetic mass has been 
calculated in the high temperature limit, where the dimensionally reduced 
model is just the 3d SU(2) pure gauge theory, and the result is [75] 

m M = O.36(2)0 2 T. (13) 

It is interesting to note that this compares well with various resummation 
methods leading to self-consistent gap equations, which at one loop yield pole 
masses in the range of rriM ~ 0.25 — 0.38g 2 T j7H], while a two-loop calculation 
gives rriM = 0.34g 2 T [77j. Note that these partonic screening masses are 
much lower than those of the singlet operators, and hence correspond to larger 
correlation lengths. This should play an important role for fragmentation 
phenomena in the plasma. 

For a constituent picture of the plasma to be appropriate, it should be 
possible to interpret the singlet operators as multigluon states For 
high T within the 3d theory, this means [75] that the masses of the lowest 
dimensional operators should be approximated by appropriate multiples of 
mg, rriM- This appears to work up to about 10% for temperatures T > 100T C 
[751 l79j. but accuracy decreases with T. E.g. the two masses shown in Figure^] 
(left) seem to slowly approach their multigluon ratio 3/2 from above |4*5] . 



6 Free Energy of Static Quarks 

A different source of information on gluonic excitations is provided by heavy 
quark free energies, in particular of a quark-antiquark pair. Such systems 
are also of interest for the physics of heavy quarkonia in the medium, in 
particular for the question of J/^-suppression jHU]. The QQ free energy is 
defined [HI] by the partition function of the thermal heat bath containing a 
static quark and antiquark source at separation R, 

( trL(R) trLt(O) ) = exp{-(F Q(5 (\R\,T) - F (T) )/T}, (14) 

where trL(x) = (l/iV c )tr Y[ T U T (x, r) is the Polyakov loop and F denotes the 
free energy of the heat bath. The gauge invariant Polyakov loop correlation 
averages over color singlet and octet contributions, and the free energy is a 
superposition [ED El EH], 

e - F QQ( R < T )/ T = - e ~Fi{R,T)/T + 8 e -F 8 (R,T)/T_ /-^ 
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Figure 8: Left: Static quark free energy for Nf = 3 at temperatures 0.58 < 
T/T c < 1.15 jHl- F(r) is normalized at (r = 1/T) to the T = Cornell 
potential, V{r)/y/a = —a/r^Ja + r^fa with a = 0.25 ± 0.05 (solid band). 
Right: The (quenched) screening mass for T > T c from fits to Eq. (J16)) with 
d = 3/2 (HZO - 



At zero temperature the free energies reduce to the heavy quark potentials, 
and at non-zero temperature they also exhibit modifications of the potentials 
by Boltzmann weighted thermal excitations. 

In the pure gluon theory, at temperatures below T c and large distances 
R, the free energy rises linearly. The coefficient of the linear term, the string 
tension, decreases with increasing temperature and vanishes above T c 
In the presence of dynamical quarks the color charges of the heavy quarks 
are screened also below T c and one observes [SSI EH ESI the expected string 
breaking, Figure |H1 The distances where the free energies become flat in R 
range from 1.5 to 1 fm, decreasing with temperature, even at (bare) quark 
masses as large as m q /T = 0.4. Note that the deviations from the zero 
temperature quenched potential set in already at distances of O(0. 5fm) for 
temperatures > 0.75T C . Normalizing the free energy to the short distance 
zero temperature potential, its large R asymptotic value rapidly decreases 
with T. 

In the deconfined phase the free energy shows the behavior of an expo- 
nentially screened potential, and is well described by 

f qq(Ri T ) _ c(T) (T)R 

T " (RT) d ' 1 J 

where c(T),d, n(T) are fit parameters. In perturbation theory, the leading 
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term originates from two-gluon exchange and predicts d = 2 and exponential 
decay with twice the Debye mass (H21- Lattice investigations jHH ESI have 
shown that this simple behavior does not apply in the temperature range 
explored. Rather, fits [Slj to Eq. (fTljj) favored d ~ 3/2 and found screening 
masses yu(T) to be compatible with the lowest color singlet 0+ + screening 
mass shown in Figures El This is not surprising: the Polyakov loop is 
a gauge invariant operator, and since it is an exponential of gauge fields 
it couples to all J PC sectors. Consequently its correlator decays with the 
lightest screening mass of the spectrum. 

In order to address Debye screening of a static quark-antiquark system, 
one needs to access the singlet free energy channel, which in the deconfined 
phase is expected to decay exponentially with the electric mass m E |H3]. Con- 
structing thermal correlators for the different channels separately in general 
requires gauge fixing. Recently it has been shown that this can be achieved 
in a manner respecting the transfer matrix, and such correlators decay with 
a gauge invariant spectrum jHTj. The octet channel was found to be repul- 
sive at short distances, as predicted by perturbation theory (EE1IE2], while 
at large distances it approaches the confining singlet potential from above. 
First studies of the thermal singlet channel have been reported in [55] . 



Susceptibilities of quark number and other quantities offer an elegant and 
simple way to study effects of small finite baryon densities on the equation of 
state within simulations performed at zero chemical potential [90J. Moreover, 
some of them are directly related to event-to-event fluctuations of particle 
production in heavy ion collisions, and help to providing signals to differen- 
tiate between the phases in order to detect the quark gluon plasma pi] . 
Quark number densities and susceptibilities are defined by 



where the jif are chemical potentials for quark number of flavor /, and rif = 
for Hf — 0. Susceptibilities for conserved charges are derived by writing down 
the appropriate superposition of quark chemical potentials and applying the 
chain rule. First we discuss the flavor singlet and triplet quark number 
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Figure 9: Left: Triplet and strange quark number susceptibilites on 12 3 x 6 
and 16 3 x 8. The lines on the right indicate free quark values for N t = 6, 8 
and the continuum [53]. Right: Response of the pseudoscalar screening mass 
to chemical potentials, 16 x 8 2 x 4, N f = 2, am = 0.025 (m n « 300 MeV) 0(i . 



susceptibilities in the two flavor case, Xsing = ^Xuu+^Xud, Xtrip = IXuu-^Xud- 
Numerical studies with dynamical staggered quarks show a rapid jump at 
the critical temperature A more recent study [93J for 2+1 flavors 

is shown in Figure El where the strange quark mass was fixed by tuning 
m riss/ m <f> to its physical value, i.e. m ! f s {2GeV) « 95 MeV [21]. In addition 
it contains the diagonal strange quark number susceptibility, which at large 
temperatures should become equal to half the triplet susceptibility. The 
use of improved actions (one loop Symanzik for gauge fields and Asqtad 
for fermions) leads to good apparent scaling behavior, suggesting deviations 
from the free quark gas value at 2T C of not more than 10 % in the continuum 
limit. 

Other studies jODl EH1 considered mostly the deconfined phase, T > 
1.5T C . They find that in this regime quenching only affects susceptibilities 
by about 5%, while cut-off effects for unimproved actions play an important 
role. Based on these observations, quenched calculations using two differ- 
ent fermion discretizations were used to obtain continuum extrapolated re- 
sults jSHl, as shown in Table 121 These authors also study susceptibilites 
for baryon number, xb = (4%3 + Xs + ^Xud + ^Xus), and electric charge, 
Xq = (10x3 + Xs + Xud ~ 2>Xus), where X3 is Xtrip /2. Flavor off-diagonal sus- 
ceptibilities are found to be consistent with zero, which is in contrast with 
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Table 2: Results for the continuum limit of quark number susceptibilities in 
quenched QCD, extrapolated from staggered quark simulations with valence 
mass m s /T c = 1. From |5^] . 



T/T c 


Xs/T 2 


Xud/T 2 


Xs/T 2 


WT 2 


xq/t 2 


1.5 


0.84 (2) 


(-2 ±3) x 10- 5 


0.53 (1) 


0.43 (1) 


0.99 (2) 


2.0 


0.89 (2) 


(-4 ±4) x 10~ 6 


0.71 (2) 


0.47 (1) 


1.07 (2) 


3.0 


0.90 (3) 


(2 ±2) x 10~ 6 


0.84 (3) 


0.49 (1) 


1.09 (3) 



predictions from a three-loop perturbative calculation [94J. It is also inconsis- 
tent with the results of [93], which reports statistically significant off-diagonal 
elements of the order of 10~ 2 at T ~ 1.5T C . While this issue needs clarifi- 
cation, off-diagonal values in any case seem to be small compared to the 
diagonal ones except very close to T c [53] . 

For the isovector susceptibility %3 numerical results are compatible. At 
T = 3T C , X3 from Table |2] falls short of its free fermion gas value Xff = T 2 
by about 10%. For such temperatures and larger, it is consistent with HTL- 
resummations, which predict Xz/Xff ~ 0-90 ~ 0.94. This observation fits 
to the behavior of the pressure and the discussion of screening masses at 
comparable temperatures, with soft gluon modes accounting for remnant 
interactions. It is interesting to note that for the isovector chemical potential 
^3 = ^11 — the fermion determinant is positive, thus it can be simulated 
(cf. Section 12]) and numerical results can be compared to analytic predictions 

In a similar vein, responses of the lowest pseudoscalar screening mass to 
scalar and vector chemical potentials have been investigated in a dynamical 
simulation with two light staggered flavors [20] • Because of the symmetry 
of the grand canonical partition function under combined euclidean time 
and \i- reflections, terms odd in (p/T) vanish for screening masses, which is 
also observed numerically. The result for the leading quadratic coefficient is 
shown in Figure 0] where [is = fi u = ^d, and \iy = \i u = —fid- Again one 
observes a pronounced jump across the deconfinement transition. Below T c , 
the response to a small scalar chemical potential is essentially zero, consistent 
with the fact that pions are still Goldstone bosons associated with chiral 
symmetry breaking in this regime. This ceases to be the case above T c , 
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where a large response is observed. The behavior with isovector chemical 
potential is quite different, and consistent with expectations based on the 
phase structure conjectured in [§o] . 

8 Dynamic Properties 

Lattice simulations inevitably are carried out in euclidean space-time. In 
order to make predictions for real-time processes, in general an analytic con- 
tinuation to Minkowski space has to be performed. At zero temperature, 
this is trivial for certain quantities such as masses of stable states or certain 
matrix elements. At finite temperature the situation is more complicated 
PP, since e.g. the momentum space propagator is defined only at discrete 
Matsubara frequencies iu n . Furthermore, because Lorentz invariance is lost 
due to the presence of the heat bath, spatial and temporal correlation func- 
tions are different in general, and the static results as discussed in Section |3] 
can not immediately be used in the interpretation of experimental results for 
dynamical processes. 

The full information about plasma excitations, the presence of genuine 
particle poles, resonances, their location and widths etc., is contained in the 
spectral density <7h{po,P) i n a given channel with quantum numbers H. It 
is related to temporal correlation functions by 



Temporal correlations have been used in order to gain information about 
the spectral density. However, since the extension of the system in the 
temporal direction is limited by the inverse temperature, the isolation of 
a ground state dominating the correlation function at large distances is dif- 
ficult. To improve the reliability of fits, extended (smeared) operators of 
various kinds have been used, sometimes combined with anisotropic lattices 
[HnilHll- With smearing one hopes to increase the projection onto the ground 
state, whose contribution would then dominate already at small temporal dis- 
tances. Anisotropic lattices increase the number of data points in the tem- 
poral direction, to stabilize more sophisticated fit ansatze including ground 
and first excited states. However, one still has to rely on fit ansatze and, 
if possible, use the quality of the fit to distinguish between various mod- 
els. Moreover, modeling the operators introduces bias, and the results have 




(18) 
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thus turned out to depend on the method. However, at least qualitatively 
the results are in agreement with a significant two-quark cut contribution, 
while in addition the behavior of wave functions was interpreted to indicate 
meta-stable bound states in the plasma phase 

Clearly, it would be much more preferable to obtain the spectral density 
u n directly. However, since the lattice provides only a limited set of noisy 
data points at discrete values of r, its numerical extraction from correlation 
functions is an ill-posed problem. Here, progress has been made recently 97 
by applying the maximum entropy method (MEM) [98J, which attempts to 
construct the most probable spectral density given the data and taking into 
account prior knowledge, such as the positivity of an and its perturbative 
behavior at large energies. The method has been successfully tested and ap- 
plied at zero temperature jHZl EH] , where the contributions of several excited 
states to the spectral density could be established. At finite temperature one 
still is hampered by the limited temporal extent of the lattice. Moreover, it 
is important to check whether one is able to reproduce the two-quark cut in 
the free case |lUUj , because this is expected to contribute significantly also in 
the interacting case at high temperature. First results for the vector channel 
at small quark masses |lUlj are shown in Figure The main feature of 
the data is that the rho peaks at low temperature disappear and develop 
into broad "resonances" above T c whose locations move proportional to the 
temperature. The vector spectral density oy(po> 0) is immediately related to 
the thermal cross section for the production of dilepton pairs at vanishing 
momentum 102 , 

dW 5a 2 1 

d^d^=° = 27^(eWT_i)^'°) • ( 19 ) 

This is shown in Figure ITU1 (right) . The "resonance" like enhancement of ay 
results in the enhancement of the dilepton rate over the perturbative tree 
level (Born) rate |lU3j for energies po/T 6 [4, 8]. In the low energy range the 
spectral density decreases rapidly, in contrast to perturbative calculations 
which have ay diverging in this limit I102L Fortunately, from inspecting 
the correlation function directly |lUlj it is already clear that the spectral 
density vanishes ~ p$ with some power a in the po ~ y limit. However, 
gaining detailed control over the low energy behavior of the spectral density 
is a challenge at this point and will become even more difficult at higher 
temperatures. On the other hand, the zero energy limit of spectral functions 
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Figure 10: Left: Vector spectral functions for various T. Right: Dilepton rate 
at 1.5 and 3T C , together with tree level and hard thermal loop prediction. 
From [TUT] . 

is related to transport coefficients |lU4j . and it would be extremely interest- 
ing to obtain results for these quantities s 105 a as a step into first principles 
investigations of non-equilibrium properties of QCD. 

9 Lattice QCD at Finite Density 

QCD at finite baryon density plays a role in two rather different regimes. 
In heavy ion collisions the initial state has a small non-zero baryon number, 
corresponding to baryon chemical potentials of order ~ O(50 MeV) |1U6] . and 
any subsequent plasma state is a state of high temperature and low density. 
On other hand, the core of neutron stars is composed of cold and very dense 
nuclear matter. These two situations correspond to the regions close to the 
axes in the tentative QCD phase diagram Figure While in the latter 
case a rich phase structure including a color superconducting state has been 
conjectured |108j . understanding the former regime is particularly pressing 
in view of current heavy ion collison experiments. Of special importance is 
a prediction of the location of the tricritical point, in which the first order 
deconfinement transition line is expected to terminate. 

Simulations of lattice QCD have so far failed to be a viable tool for finite 
density because of the so-called "sign problem". For \i = the relation 
75M75 = M< guarantees positivity of the fermion determinant, det M(0) > 0. 
However, for SU(3) and /1 / this relation does not hold and detM is 
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Figure 11: Conjectured phase diagram for QCD with two massless and one 
heavy flavour. From |107j . 



complex. This prohibits Monte Carlo importance sampling, which interprets 
the measure as a probability factor and thus requires it to be positive. Note 
that all physical quantities like Z or observables (O) are real positive, i.e. all 
imaginary parts cancel out of the path integral and the problem lies merely 
in evaluating it. 

This problem is generic for many fermionic models, and for a certain class 
of them it can actually be solved by the use of cluster algorithms. These 
subdivide the ensemble into parts within which all contributions add up 
positively. The model closest to QCD for which this succeeded is the three 
state Potts model |109j . Unfortunately, these algorithms are not directly 
applicable to QCD, and we refer to |110j for an overview. 

We will likewise not discuss the considerable amount of work done on 
certain modifications of QCD for which the sign problem is manageable or 
altogether absent, such as QCD in the static quark limit, two-color QCD or 
QCD at finite isospin. These topics are covered in recent reviews |lllj . It 
has also been argued that for very high densities an effective theory free of 
the sign problem can be derived j!12j . but no numerical results exist as of 
yet. 

Here we wish to focus on some recent approaches which made significant 
progress towards simulating QCD with realistic parameter values at small 
baryon densities. It is important to note that all of them circumvent the 
sign problem rather than solving it. The partition function in its first form 
in Eq. (J2J suggests that for small £l = [i/T the problem is close to the trivial 
case, and in particular lends itself to a Taylor expansion. The following 
approaches in one way or another all require small enough Ji to succeed. 
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As we shall see, however, this limitation still allows to obtain results up to 
fiB ~ 500 MeV, thus fully including the parameter region relevant for heavy 
ion collisions. In the following we limit our discussion to mapping out the 
deconfmement transition line. 

9.1 Reweighting: The Glasgow Method 

A mathematical identity can be used to reweight from configurations gener- 
ated at zero density, where there is no sign problem, to a corresponding one 
at finite density. The method is based on rewriting the partition function as 

Z = fpU detM(O) P MM e -».W = , (20) 

J v; detM(0) \det(M(0))/ ' V ' 

where the determinant ratio is now treated like an observable while the in- 
tegration measure is defined at /i = 0, and hence positive. This expression 
is exact. However, in a Monte Carlo integration with finite statistics two 
problems arise. Firstly, the complexity of the numerator leads to oscilla- 
tions and large cancellations. The reweighting factor corresponds to a ratio 
of two partition functions with different actions, and thus decays exponen- 
tially with the difference between their free energies, which is an extensive 
quantity, Z = exp — AF/T ~ exp (—const. V). The statistics required for a 
given accuracy thus grows exponentially with volume, rendering extrapola- 
tions to the thermodynamic limit extremely difficult. Secondly, Monte Carlo 
approximates the partition function by the contributions of the most likely 
configurations obtained by probability sampling. Since the extrema of the 
action with /i = are shifted from those of the entire integrand including the 
determinant ratio, the probability distribution of the generated ensemble is 
changed. This should be of little consequence for small /i and large statis- 
tics, but with growing /i the overlap between reweighted and full ensemble 
deteriorates. The effect does not show up in the statistical error, and the 
difficutly is in knowing when this becomes problematic leading to incorrect 
results. For example, at T = the transition to nuclear matter is found 
at an unphysical value for /i c . A detailed discussion with early results and 
references can be found in |113| . 
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9.2 Multiparameter Reweighting 



Only recently it was realized that the ensemble overlap can be significantly 
improved by a multidimensional generalization |114j of the Glasgow method. 
In addition to /i one may also reweight in the lattice gauge coupling (3, so 
that the partition function is now written as 



This is crucial for finding critical behavior. While reweighting only in /i would 
inevitably generate an ensemble away from criticality, the second reweight- 
ing parameter can be used to keep it fluctuating between the phases, as the 
unshifted ensemble certainly would. This method was compared with the 
Glasgow method in four flavor QCD at imaginary chemical potential |114j . 
for which the determinant is positive and a full ensemble is also available 
(cf. Section inHJ). Multiparameter reweighting is found to agree significantly 
better with the full ensemble up to some critical imaginary chemical poten- 
tial, beyond which it breaks down |115j . 

The method was then used in to map out the deconfinement line, 
using Lee Yang zeros |117j as observables. These are points of vanishing par- 
tition function, marking the singularities in the free energy. Singular phase 
transitions only appear in the thermodynamic limit, while on finite volumes 
the free energy is regular. In the latter case the Lee Yang zeroes are displaced 
to complex parameter values. As the volume is increased, they approach their 
thermodynamic limit values on the real axis if there is a genuine phase tran- 
sition, whereas for a smooth crossover the zeroes settle somewhere in the 
complex plane. In the case of a first order transition the large volume be- 
havior is consistent with a (3 C (V) = /3 c (oo) + (/V behavior. The simulations 
reported in |116j were carried out for 2+1 flavors of staggered fermions, the 
result is shown in Figure It constitutes the first numerical prediction of 
the {n,T) phase diagram including a critical point. However, larger volumes 
are required to check for thermodynamical behavior and confirm the results. 
Recently attempts have been made to further improve the overlap by split- 
ting the reweighting factor in a product of many factors, and preliminary 
results on small lattices have been reported in |118j . 




Z 



e -s g (P) det(M(/x)) 
e -5 9 (A,)det(M(0)) 



(21) 
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9.3 Taylor Expanded Reweighting 

Since for heavy ion collisions one is interested in rather small chemical poten- 
tials of a few ten MeV, one may use the smallness of ft in an approximation of 
multiparameter reweighting in order to make simulations on larger volumes 
feasible |119j . Consider a Taylor expansion of the reweighting factor TZ as 
a power series in ft = [i/T, and similarly for any operator O. Expectation 
values are then given by 

{ )iM (ex V (n 1 ft + n 2 ft 2 + ...-As g )) p=0A 

(22) 

For sufficiently small chemical potentials, this series should converge quickly 
and calculation of the leading coefficients should give a good approximation 
to the full answer. The benefit is that the derivatives of the determinant 
can be expressed through local operators in the fermion fields. These are 
much cheaper to compute than the whole determinant, which requires a 
non-local operation. In particular, keeping only the leading term, the cost is 
equivalent to that of computing susceptibilities (cf. Section Ej). Even though 
the expanded phase of the determinant still has oscillations growing with 
volume and ft, the saving in computer time can be used to gather higher 
statistics for larger volumes. 

This strategy has been applied in |119j to compute the location of the de- 
confinement line for a two-flavor theory with p4-improved staggered fermions 
and a pion mass of about « 600 MeV. The method to locate a phase 
transition is to look for the peak Xmax = xif^c, Pc) of the susceptibilities 
X = VN t ((O — (O)) 2 ), where the chiral condensate and the Polyakov loop 
were used for O in practice. This peak implicitly defines a critical coupling 
Pdft), which must be even in \i because of the symmetry of the partition func- 
tion This can be employed to calculate the curvature of the critical 
temperature at /i = 0, 



d 2 T c 1 d 2 (3 c 



dn 2 N 2 T C d\i 2 




(23) 



Indeed a lattice with twice the spatial size than that in was employed 
without major problems. Of course, in order to judge the convergence prop- 
erties of the series, several consecutive coefficients are needed. 
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9.4 Imaginary Chemical Potential 



For imaginary chemical potential, the fermion determinant is positive and 
simulations are as straightforward as for \x = 0. It is therefore natural to ask 
whether such simulations can be exploited to learn something about real \i. 
Let hr, Hi e 1R denote the real and imaginary parts of a complex \i = hr+' 1 Hi- 
In the presence of a complex chemical potential, a Z(3) transformation of 
the fermion fields is equivalent to a shift in hi, causing the partition function 
to be periodic in Hi with period 2ttT/3 |120j . 



Hence, once Hi exceeds some critical value flh a phase transition to a non- 
trivial Z{3) sector occurs, and because of the symmetry Eq. fl21| this transi- 
tion is periodically repeated at the exact critical values fij = 2n(k + l/2)/3. 
It was conjectured |12Uj and verified numerically |115H12lj that these transi- 
tions are of first order in the deconfined phase and continuous in the confined 
phase. On the other hand, for purely real \i these transitions are absent. 

Through its periodicity, the grand canonical partition function at imag- 
inary chemical potential is the Fourier transform of the canonical partition 
function at fixed quark number Q |12U| I122j , 



It was thus suggested to simulate Z at imaginary chemical potential and 
perform the Fourier transformation numerically |123| . This has been car- 
ried out in the 2d Hubbard model at large T and small Q, but not yet in 
QCD. Unfortunately, this approach only postpones the sign problem from 
the Monte Carlo integration to the Fourier transform. As Q is getting large 
approaching the thermodynamic limit at fixed baryon density, it will again 
cause rapid oscillations and the associated cancellations. 

9.5 Analytic Continuation 

If fi is small another line of attack is possible |124j . Since there are no massless 
modes in the theroy, all observables are analytic functions of fi everywhere 
but on the critical lines of phase transitions. After computing expectation 
values with Z(i/xj,T, V), one may therefore fit them by truncated Taylor 



Z(hr, pi) = Z(jiR, Hi + 27rT/3). 



(24) 



Z Q (T, V) = — dni e-'^ Q ' T Z{H = m, T, V) 

ZTT JO 



(25) 
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series 



A' 



(26) 



n 



In cases where this is possible to satisfactory accuracy, analytic continuation 
to real /x is trivial. There are two advantages to such a procedure. First, 
there is no sign problem and hence no restriction on the volumes that can be 
simulated. Second, the simulation gives the value of the whole expectation 
value rather than just a single Taylor coefficient. While the expansion is 
necessary in order to analytically continue, fitting to various orders allows 
control over the convergence properties. However, a check of the convergence 
and hence controlled continuation are limited to /2j < p,j — 7r/3, which marks 
the first Z(3) transition discussed above. 

Non-perturbative evidence for the viability of this approach has been 
given for screening masses in the deconfined phase |125| . which can be simu- 
lated sucessfully for imaginary and real /i in the framework of dimensionally 
reduced QCD |52|. 

It was shown in |115j that the same approach can be extended to study 
the critical line of the deconfinement transition itself. Again the transition 
was defined by the peak Xmax = xif^a A) °f plaquette, chiral condensate or 
Polyakov loop susceptibilities. On finite volumes these are analytic functions 
over the whole parameter space of the theory, non-analyticities associated 
with phase transitions develop only in the thermodynamic limit. The implicit 
function theorem then implies that also A Of) is analytic for all p,, which 
therefore has a Taylor expansion 



As the thermodynamic limit is approached, A(aO approaches the unique 
critical line of phase transitions, while remaining pseudo-critical in a crossover 
regime. In principle the nature of the line can also be determined from its 
finite volume scaling. 

In |115j a two staggered flavor QCD was simulated, and a striking finding 
is that for all /jj < p,} = 7r/3 the critical coupling is indeed well described 
by a fit quadratic in p,j, with negligible effect of terms O(flj). The same 
conclusion is reported from an imaginary \x simulation of a four flavor theory 
in [121] . In principle, of course, there is nothing that guarantees such good 
behavior, which would need to be checked for every change of parameters. 
It is therefore amusing to note that both in the 3d Gross-Neveu model as 



(27) 



n 
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well as in some random matrix model with symmetries similar to QCD, the 
leading quadratic term very well reproduces the exact soultion for the critical 
line over wide ranges of parameter space |121j . 

9.6 The (/i, T) Phase Diagram 

In Figure ED we assemble the various phase diagrams published so far in 
one plot. Some caution shoud be used when interpreting it, as the curves 
correspond to different parameter values and volumes. The only thing equal 
in all four simulations is the lattice spacing which at a ~ 0.3 fm is large 
enough to cause cut-off effects. Furthermore, the error bars shown on the 
data points should not be directly compared with the smaller error bands 
corresponding to fit results. 




200 400 600 

|i /MeV 



Figure 12: Location of the (pseudo-) critical line as reported in four different 
simulations. All employ staggered fermions, only Allton et al. use improved 
actions. The simulation parameters are collected in Table El 

Nevertheless, the curves are only shown as far as the validity of the 
quadratic approximation has been checked explicitly in the imaginary /x ap- 
proach. Within this range it was found in |119j that doubling the quark 
mass did not affect the line within the current errors, so that comparing 
different quark mass values makes sense. Moreover, the four flavor calcu- 
lation at imaginary fi in |121j is entirely consistent with the (unexpanded) 
reweighting calculation performed at exactly the same parameter values |114j . 
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Table 3: Simulation parameters for the results shown in Figure IT21 



Method 


N f 


m q 


largest lattice 


reweighting |116j 


2+1 


am = 0.025, m s = 8m u , {m* w 300Mel/) 


8 3 x 4 


rew. + Taylor |119j 


2 


am = 0.1, (m* « 600) MeV 


16 3 x 4 


imag. fi [115J 


2 


am = 0.025, (m* « 300) MeV 


8 3 x 4 


imag. jx j!21j 


4 


am = 0.05 


16 3 x 4 



Finally, within the Taylor expanded reweighting it was checked that simulat- 
ing real and imaginary chemical potential indeed produces the same leading 
coefficient [126J. 

These cross-checks are extremely important, since all methods introduce 
some systematic error growing with fi/T. However, it seems safe to conclude 
that they all give compatible and reliable results for the pseudo-critical tem- 
perature at least up to \i& ~ 500 MeV. As expected on physical grounds, 
T c (/i) decreases faster when more light quarks are present. An important 
physics conclusion then is the flatness of the critical line, which implies that 
the critical temperature relevant for RHIC is essentially the same as at zero 
density. It also implies that T c (fis —* lGeV) has to drop rather dramatically 
if our qualitative picture of the phase diagram is at all correct. First inves- 
tigations of the equation of state at finite /i have found negligible effects at 
the few percent level for RHIC densities |127| I119j . 

9.7 The Critical Point, Quark Mass Dependence 

The only simulation that has so far produced a critical point is the one em- 
ploying reweighting on lattices of rather limited size |116j . On the other 
hand, experience at zero density indicates that rather large volumes are re- 
quired in order to reach the asymptotic finite volume scaling required for 
an unambiguous determination of the order of a phase transition |19j . Such 
systematic effects may well be larger than the statistical error box indicated 
in the figure, and it is hence important to have another independent calcu- 
lation. Moreover, while the location of the (pseudo-) critical line appears to 
be only very weakly quark mass dependent, the location of the critical point 
strongly depends on them: for three chiral quarks, the first order line joins 
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the temperature axis, while for the 2+1 flavor case depicted in the figure it 
is quite far out. 

The task of locating the critical point as a function of parameters is 
equivalent to mapping out the critical surface separating first order transi- 
tions from crossover in Figure which is simply the extension of the phase 
diagram Figure Q to finite densities. Probably the most economical way of 
doing this is to find the projections of this surface onto the (T, /i) and various 
(m q , //)-planes. Just like the former, the (m q , yu)-plane then shows a critical 
line separating first order from crossover. For this line, m q (fi), the same an- 
alyticity and symmetry considerations apply as discussed for T c (/z), i.e. it 
has a Taylor series even in /i. Within the Taylor expanded reweighting ap- 
proach, a first result for the leading coefficient in the case of three degenerate 
flavors of p4-improved fermions (m^ w 200MeV) has been reported in [23] 
as Tdm c q / d(ii 2 ) = 0.21(6). Moreover, as argued in |128j . when the quark 
mass is varied, the 3 flavor critical endpoint will trace out a line T c (/x c , m q ), 
which has an analytic continuation to imaginary /i. This in turn allows to 
determine it from imaginary chemical potential simulations. 




Figure 13: Critical surface separating the QCD parameter space for first 
order transitions from crossover. 



10 Summary 

In this article we have attempted to summarize the current understanding 
of QCD at non- vanishing temperature as it arises from numerical studies on 
the lattice. While many of the results in the pure gauge theory are available 
in the continuum limit, simulations with dynamical fermions still suffer from 
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systematic errors. These are mainly due to finite lattice spacing as well as 
quark masses which don't meet their physical values yet. Nevertheless, for 
critical temperatures and the equation of state quantitative results could be 
given, which are expected to be correct at the 10 - 15 % level. 

Apart from the basic thermodynamic quantities, existing results provide 
us with a detailed picture of how static quantities change through the decon- 
finement transition up to a few T c . Combinations of perturbative calculations 
and numerical simulations have produced insight into the regime of very high 
temperatures as well as the dynamics and mixing of modes. Altogether this 
led to a quantitative understanding of the relevant static length scales in 
the plasma, as well as tests of the applicability and breakdown of thermal 
perturbation theory. 

The naive picture of the deconfined phase as a weakly interacting par- 
ton gas is not supported. For temperatures relevant to heavy ion collisions, 
the plasma displays strong residual interactions through soft gluonic modes, 
which cannot be treated perturbatively, and which influence different quanti- 
ties in different ways. In particular, any constituents themselves are objects 
dressed by non-perturbative interactions. This gives a consistent explana- 
tion to the various observed features: the equation of state, susceptibilities 
and fermionic correlators are dominated by hard modes, but significant de- 
viations from ideal gas behavior are still present. Gluonic correlators, on 
the other hand, are dominated by soft modes and entirely off their leading 
perturbative predictions. An ideal gas is established only at asymptotically 
high temperatures. 

The study of real time dynamics in the plasma is a much more difficult 
problem, for which at present no formal approach exists within lattice gauge 
theory. A step towards inspecting such physics could be presented by means 
of maximum entropy methods. However, at this stage such calculations have 
still a qualitative character, in trying to establish a working methodology. 

A large section of this paper was devoted to another new development, 
allowing to study QCD also at small baryon density, as it is explored by heavy 
ion experiments at RHIC and LHC. Again, the presented studies are still at 
an exploratory level, subject to the already mentioned systematic effects. 
However, rather different methods have led to a consistent dependence of 
the critical temperature on the chemical potential, giving confidence in their 
validity up to /ib ~ 500 MeV. In this range the pseudo-critical line is found 
to be rather flat. There also is a prediction for a critical point to be cross 
checked in simulations to come. 
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Advances in lattice QCD have always been brought about by a mixture of 
new calculational schemes, new algorithms and refined analysis techniques. 
They have, however, also been related to progress in computational resources. 
With a new generation of machines in the Teraflops-range going into opera- 
tion soon, one has all reasons to hope that part of the mentioned systematic 
deficiencies will be reduced and simulations at parameter values much closer 
to the physical ones will be realized. We hope to have also shown that, by 
combining analytic and numerical techniques, insight can be gained and new 
paths may show up, which might be needed to venture fully into the real 
time regime. 
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